/* **************************************************** */
/* * Gauss Method for finding PI - Sequential Version * */
/* **************************************************** */

#include <gmp.h>
#include <stdio.h>
#include <stdlib.h>

/* Number of iterations */
#define NUM_ITERS 2000
/* Size of array. For best results, NUM_ITERS % ARRAY_SIZE should be 0 */
#define ARRAY_SIZE 2000

/* Global array of mpf_t variables */
mpf_t a[ARRAY_SIZE+1], b[ARRAY_SIZE+1], p[ARRAY_SIZE+1], tn, tnp;

int main(int argc, char *argv[]) {
	int k, i, errorNumber = 0;
	mpf_t pi, piBottom;

	/* Define precision */
	if(2 == argc)
		mpf_set_default_prec(atoi(argv[1]));
	else if(argc > 2) {
		printf("Usage: %s [prec]\n", argv[0]);
		exit(1);
	}

	/* Set a[0] = 1 */
	mpf_init_set_d(a[0], 1.0);
	/* Set b[0] = 1/sqrt(2) */
	mpf_init_set_d(b[0], 2.0);
	mpf_sqrt(b[0], b[0]);
	mpf_ui_div(b[0], 1, b[0]);
	/* Set p[0] = 1 */
	mpf_init_set_d(p[0], 1.0);
	/* Set t[0] = 1/4 */
	mpf_init_set_d(tn, 0.25);
	/* Initialise t[i+1] */
	mpf_init(tnp);
	/* Initialise PI values */
	mpf_init(pi);
	mpf_init(piBottom);

	/* Initialise the rest */
	for(i = 1; i <= ARRAY_SIZE; i++) {
		mpf_init(a[i]);
		mpf_init(b[i]);
		mpf_init(p[i]);
	}

	/* Start process */
	for(k = 0; k < NUM_ITERS; k += ARRAY_SIZE) {
		for(i = 0; i < ARRAY_SIZE; i++) {
			/* Calculate a[i+1] */
			mpf_add(a[i+1], a[i], b[i]);
			mpf_div_ui(a[i+1], a[i+1], 2);

			/* Calculate b[i+1] */
			mpf_mul(b[i+1], a[i], b[i]);
			mpf_sqrt(b[i+1], b[i+1]);

			/* Calculate p[i+1] */
			mpf_mul_ui(p[i+1], p[i], 2);

			/* Calculate t[i+1] */
			mpf_sub(tnp, a[i], a[i+1]);
			mpf_pow_ui(tnp, tnp, 2);
			mpf_mul(tnp, tnp, p[i]);
			mpf_sub(tn, tn, tnp);
		}

		/* Recycle array: Last elements become first and process is restarted */
		mpf_set(a[0], a[ARRAY_SIZE]);
		mpf_set(b[0], b[ARRAY_SIZE]);
		mpf_set(p[0], p[ARRAY_SIZE]);
	}

	/* Calculate PI using final values */
	mpf_add(pi, a[0], b[0]);
	mpf_pow_ui(pi, pi, 2);
	mpf_mul_ui(piBottom, tn, 4);
	mpf_div(pi, pi, piBottom);

	/* Print final results */
	gmp_printf("Final a[n]: %.6Ff\n", a[0]);
	gmp_printf("Final b[n]: %.6Ff\n", b[0]);
	gmp_printf("Final p[n]: %.6Ff\n", p[0]);
	gmp_printf("Final t[n]: %.6Ff\n", tn);
	gmp_printf("Final PI: %.6Ff\n", pi);

	/* Free everything */
	for(i = 0; i <= ARRAY_SIZE; i++) {
		mpf_clear(a[i]);
		mpf_clear(b[i]);
		mpf_clear(p[i]);
	}
	mpf_clear(tn);
	mpf_clear(tnp);
	mpf_clear(pi);
	mpf_clear(piBottom);

	return 0;
}